\documentclass{ctexart}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{subfigure}
\usepackage{array}
\usepackage{float}

\title{微分方程数值解hw02 \\ Chapter7-Programming Assignments}


\author{陈冬勇 \\ 信息与计算科学 \ 3200103953}

\begin{document}

\maketitle

\section{编译说明}

\noindent
本次项目作业由\ \verb|Makefile|文件进行编译，以\ \verb|FDmethod.h|为头文件，输入\ \verb|make run|便能由\ \verb|main.cpp|产生可执行文件\ \verb|test|（一些数值、函数输入都已经存放在程序内部以及.json文件当中）并同时生成tex的.pdf文件，输入\ \verb|make clean|便能清理多余的文件，输入\ \verb|./test|便能运行可执行文件并得到结果存入对应的.txt文件，利用这些文本文档可以在python中作图得到.png文件用于报告内。实验报告中需要引用课本中定义、定理、引理、式子等会使用粗体字体进行表示。 \\

\section{FDmethod.h头文件说明}

\noindent
在头文件中，一共设计了Function和FD两个基类，派生出不同的class以满足项目要求。（鉴于本人能力有限，解决不了第二种区域的问题，有关FD的class FD\_domain2以及其所需的有关圆的类都没有写在头文件内） \\
Function类：一共设计了operator()，diff\_x，diff\_y和Laplace四个虚函数，用于返回给定二元函数的值、一阶偏导数的值以及负Laplace算子作用后的值。其主要功能是存储待求解函数u的表达式，以便输入边界条件并对其进行一系列的误差分析。 \\
FD类：FD类一共设计了solve，result和norm\_error三个虚函数，分别用于进行FD算法，返回计算得到的格点处的数值解以及返回误差的1、2以及无穷范数。其主要功能是实现基于FD方法的微分方程数值解的求解；通过继承这个类，可以创建针对单位正方形上的偏微分方程的具体算法。 \\
(1)class FD\_domain1:支持单位正方形上微分方程在Dirichelet、Neumann和Mixed conditions三种边界条件下求解。其可以接收函数值、网格宽度和边界条件等参数。solve函数求解AU=F的线性方程组，其中A为系数矩阵并且F为向量。我们加入边界，并且按x轴正方向、y轴正方向的顺序给格点标注序号以便后续计算。该函数在满足不同的边界条件下，使用不同的方法来构造A和F: 当condition=1时，表示Dirichlet边界条件(简称D条件)。此时，在矩阵A的内部区域中，使用五点差分公式来近似求解二阶偏微分方程。在边界上，将u的值作为已知量，构造出系数矩阵A和向量F；当condition=2时，表示Neumann边界条件(简称N条件)。此时，在矩阵A的内部区域中，仍然使用五点差分公式来近似求解二阶偏微分方程。而在边界上，我们根据FD算法得到的二阶精度的导数来作为导数值（左右边界为关于x的偏导数，上下边界为关于y的偏导数）的估计；当condition=3时，表示Mixed边界条件(简称M条件)，在本程序中假设与x轴平行的两边为D条件，另外两条边为N条件，此时结合以上两种方法即可构建方程组。 \\
(2)class FD\_domain2:未解决，仍待后续研究…… \\

\section{输入输出说明}

\noindent
(1)main.json输入参数说明（参数可灵活调整）： \\
n：单位正方形格点等分数。根据题目要求n=8,16,32,64。\\
function\_number：函数标号。共有三个函数需进行求解等操作，1表示给定函数$u=e^{y+sin(x)}$，2表示自定义函数$u=sin(xy)$，3表示自定义函数$e^{−xy}$。 \\
boundary\_number：边界条件标号。1表示D条件，2表示N条件, 3
表示M条件。 \\
(2)main.cpp输出结果说明： \\
只需运行可执行文件test，即可得到一个.txt文件，其中记录了格点函数值、各种范数的误差等数值用于后续的绘图和数据分析。 \\

\section{数据结果分析}

\noindent
(1)函数求解：result可以得到每一个格点处的u的数值解并将其输出至results.txt中，据此可以绘制出相应的图像（正方形区域之上的三维图像）；以下会节选某一函数某一边界条件在所有n下的图像进行呈现，共给出三种情况。 \\
A.我们取定function\_number=1以及boundary\_number=1的情况。\\
\begin{figure}[ht!]
    \centering
    \subfigure[n=8]{
    \includegraphics[width=2in]{1.png}}
    \quad
    \subfigure[n=16]{
    \includegraphics[width=2in]{2.png}}
    \quad
    \subfigure[n=32]{
    \includegraphics[width=2in]{3.png}}
    \quad
    \subfigure[n=64]{
    \includegraphics[width=2in]{4.png}}
    \caption{$u=e^{y+sin(x)}$,D条件}
\end{figure}\\
B.我们取定function\_number=2以及boundary\_number=3的情况。\\
\begin{figure}[ht!]
    \centering
    \subfigure[n=8]{
    \includegraphics[width=2in]{5.png}}
    \quad
    \subfigure[n=16]{
    \includegraphics[width=2in]{6.png}}
    \quad
    \subfigure[n=32]{
    \includegraphics[width=2in]{7.png}}
    \quad
    \subfigure[n=64]{
    \includegraphics[width=2in]{8.png}}
    \caption{$u=sin(xy)$,N条件}
\end{figure}\\
C.我们取定function\_number=3以及boundary\_number=2的情况。\\
\begin{figure}[ht!]
    \centering
    \subfigure[n=8]{
    \includegraphics[width=2in]{9.png}}
    \quad
    \subfigure[n=16]{
    \includegraphics[width=2in]{10.png}}
    \quad
    \subfigure[n=32]{
    \includegraphics[width=2in]{11.png}}
    \quad
    \subfigure[n=64]{
    \includegraphics[width=2in]{12.png}}
    \caption{$e^{−xy}$,M条件}
\end{figure}\\
(2)固定点收敛速率：在程序运行时记录求解结果与真实解在四个点 ($\frac{1}{8}, \frac{1}{8}$)、($\frac{1}{8}, \frac{7}{8}$)、
($\frac{7}{8}, \frac{1}{8}$) 和 ($\frac{7}{8}, \frac{7}{8}$) 处的绝对误差并输出至results.txt中；为分析收敛速率，首先在main.json中n的取值内加入128，增加网格数以便观察趋势；其次，求出绝对误差序列前一项与后一项的比值并取其以2为底的对数，以便观察n翻倍后绝对误差缩减的趋势。对于三个函数，我们分别记录每个函数对应的三种
边界条件的绝对误差序列的比值。（n的五个取值对应四个绝对误差比值）\\
\begin{figure}[ht!]
    \centering
    \subfigure[$u=e^{y+sin(x)}$]{
    \includegraphics[width=2in]{13.png}}
    \quad
    \subfigure[$u=sin(xy)$]{
    \includegraphics[width=2in]{14.png}}
    \quad
    \subfigure[$e^{−xy}$]{
    \includegraphics[width=2in]{15.png}}
    \caption{domain1，绝对误差比值}
\end{figure}\\
由图可知，三种边界条件下各个函数的收敛速度关于h都是2阶的。总的来说，D条件和M条件下的绝对误差比值更接近2，更为稳定；而N条件下的收敛阶并不稳定，但仍然逐渐趋于2。除此之外，根据results.txt的输出，D条件和M条件的绝对误差本身在量级上也更小，这可能与Neumann边值问题的不适定性有关。本程序针对这种特殊情况新增了一个条件以确保该偏微分方程有唯一解，不过该条件并不强，这使得数值解的误差仍然偏大。 \\
(3)误差范数及其收敛阶分析：results.txt中还会输出不同函数、边值条件和n下对应误差的1、2和无穷范数，参照$\mathbf{Definition 7.13}$即可。（注意：此处要将课本中的$h$改成$h^2$）查看.txt文件后，可以看出在多种情况下，误差的三类范数都是较小的数值（小于$1 \times 10^{-3}$），一定程度上佐证了该算法的正确性。而关于误差范数的收敛分析，首先可以理论上验证，如下所示 \\
$\Vert g \Vert_\infty=O(h^2)$ \\
$\Vert g \Vert_1=h^2(O(1))=O(h^2)$ \\
$\Vert g \Vert_2=[h^2(O(h^2))]^\frac{1}{2}=O(h^2)$ \\
其次我们可以根据.txt中的结果作图进行实验上验证，首先仍然在main.json中n的取值内加入128，增加网格数以便观察趋势；其次，求出误差范数序列前一项与后一项的比值并取其以2为底的对数，以便观察n翻倍后误差范数缩减的趋势。对于三个函数，我们分别记录每个函数对应的三种边界条件的误差范数序列的比值（n的五个取值对应四个误差范数比值），并取最简单的D条件为例进行作图。 \\
\begin{figure}[ht!]
    \centering
    \subfigure[$u=e^{y+sin(x)}$]{
    \includegraphics[width=2in]{16.png}}
    \quad
    \subfigure[$u=sin(xy)$]{
    \includegraphics[width=2in]{17.png}}
    \quad
    \subfigure[$e^{−xy}$]{
    \includegraphics[width=2in]{18.png}}
    \caption{domain1，误差范数比值，D条件}
\end{figure}\\
由图可知，误差的三种范数的收敛阶都稳定趋向于2，于是通过计算机在实验上进行了验证，并且对N条件、M条件作图也能得到类似的结果。 \\

\end{document}
